High Resolution, Arbitrary-Even-Order Direction Finding Method and Device

ABSTRACT

Method of high-resolution direction finding to an arbitrary even order, 2q (q&gt;2), for an array comprising N narrowband antennas each receiving the contribution from P sources characterized in that the algebraic properties of a matrix of cumulants of order 2q, C 2q,x (l), whose coefficients are the circular cumulants of order 2q, Cum[x i     1   (t), . . . , x i     q   (t), x i     q+1   (t)*, . . . , x i     2q   (t)*], of the observations received on each antenna, for cumulant rankings indexed by l, are utilized to define a signal subspace and a noise subspace.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present Application is based on International Application No.PCT/EP2006/061252, filed on Mar. 31, 2006, which in turn corresponds toFrench Application No. 05 03180, filed on Apr. 1, 2005, and priority ishereby claimed under 35 USC §119 based on these applications. Each ofthese applications are hereby incorporated by reference in theirentirety into the present application.

FIELD OF THE INVENTION

The present invention relates to a method of high-resolution directionfinding to an arbitrary even order, that is to say exclusively utilizingthe information contained in the statistics of order 2q of theobservations, where q is an integer such that q≧2, preferably q>2, forspace, pattern and/or polarization diversity arrays.

BACKGROUND OF THE INVENTION

The invention finds its application notably in all systems which requirea significant degree of spatial resolution [(environments with a largenumber of transmitters, array with a low number of sensors, significantresidual model errors (standardization, calibration, etc.), etc.)].

At the start of the 80s, numerous order 2 direction finding procedurestermed high-resolution (HR) [1] [14] were developed to alleviate thelimitations of procedures that were termed super resolution [2-3] inrelation to weak sources. Among these HR procedures, the so-calledsubspace procedures such as the MUSIC (or MUSIC-2) procedure [14] arethe most popular. These HR procedures are very efficacious in amulti-source context insofar as they possess, in the absence of modelerrors and for background noise of known spatial coherence, unlimitedasymptotic separating ability whatever the signal-to-noise ratio (SNR)of the sources. However, these HR procedures suffer from seriousdrawbacks. Specifically, they can process at most N−1 noncoherentsources on the basis of an array with N antennas and are not very robusteither to model errors [10] [15], which are inherent in operationalimplementations, or to the presence of background noise of unknownspatial coherence [12], typical of the HF range for example.Furthermore, their performance can become greatly affected when thisinvolves separating several weak and angularly close sources on thebasis of a limited number of observed samples.

From the end of the 80s, mainly to alleviate the previous limitations,order 4 high-resolution direction finding procedures [13] have beendeveloped for non-Gaussian sources, that are omnipresent inradiocommunications, among which the extension of the MUSIC procedure toorder 4 [13], called MUSIC-4, is the most popular. Specifically, theorder 4 procedures are asymptotically robust to the presence of Gaussiannoise of unknown spatial coherence [13]. Furthermore, in spite of theirgreater variance, they generate a virtual increase in the aperture ofthe array and the number of antennas, introducing the notion of virtualarray (RV) to order 4 [4][6] and offering increased resolution and thepossibility of processing a greater number of sources than the number ofantennas. In particular, on the basis of an array with N antennas, theMUSIC-4 procedure can process up to N(N−1)+1 sources when the antennasare identical and up to (N+1)(N−1) sources for different antennas.

The invention relates to a method of high-resolution direction findingto an arbitrary even order, 2q with q>2, for an array comprising Nnarrowband antennas each receiving the contribution from P sourcescharacterized in that the algebraic properties of a matrix of cumulantsof order 2q, C_(2q,x)(l), whose coefficients are the circular cumulantsof order 2q, Cum[x_(i) ₁ (t), . . . , x_(i) _(q) (t), x_(i) _(q+1) (t)*,. . . , x_(i) _(2q) (t)*], of the observations received on each antenna(for cumulant rankings indexed by l) are utilized to define a signalsubspace and a noise subspace.

This method presents notably the advantage of increasing virtually, bycreating virtual antennas, the aperture of the antenna arrays used andtherefore of increasing the resolution, the processing capacity and therobustness to model errors of direction finding techniques in a mannerthat increases with q.

For rectilinear sources, the method according to the invention allowsthe creation of virtual antennas and aperture on the basis of the matrixof all the cumulants of order 2q of the observations, this allowing theprocessing of a number of sources that increases with q, on the basis ofa given array of antennas.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention is illustrated by way of example, and not bylimitation, in the figures of the accompanying drawings, whereinelements having the same reference numeral designations represent likeelements throughout and wherein:

FIG. 1 the representation of an incident signal in three dimensions, and

FIG. 2 a virtual array of order 4 of a uniform circular array of 5antennas with the order of multiplicity of the virtual antennas.

DETAILED DESCRIPTION OF THE INVENTION

Before detailing the steps implemented by the method, some data usefulto an understanding of the invention is recalled.

We consider an array of N narrowband (NB) antennas, that are assumedidentical initially, and we denote by x(t) the vector of the complexenvelopes of the signals output by the antennas. Each antenna receivesthe contribution from P stationary centered NB sources, statisticallyindependent or otherwise, and background noise. It is assumed that the Psources can be grouped into G sets, with P_(g) sources in the set g,such that within a group, the sources are assumed to be statisticallydependent, while the sources belonging to different sets are assumed tobe statistically independent. In particular, G=P corresponds to the casewhere the P sources are statistically independent while G=1 correspondsto the case where all the sources are correlated. Of course, theparameters P_(g) are such that:

$\begin{matrix}{{{P = \sum\limits_{g = 1}^{G}};};P_{g}} & (1)\end{matrix}$

Under these assumptions, the observation vector can be writtenapproximately

$\begin{matrix}\begin{matrix}{{{{x(t)} \approx \sum\limits_{i = 1}^{P}};};{{{{m_{i}(t)}{a\left( {\theta_{i},\phi_{i}} \right)}} + {v(t)}} =};{{\Delta \; {{Am}(t)}} + {v(t)}}} \\{{{= \sum\limits_{g = 1}^{G}};};{{A_{g}{m_{g}(t)}} + {v(t)}}} \\{{{= \sum\limits_{g = 1}^{G}};};{{x_{g}(t)} + {v(t)}}}\end{matrix} & (2)\end{matrix}$

where ν(t) is the noise vector, assumed centered and Gaussian, m(t) isthe vector whose components m_(i)(t) are the complex envelopes of thesources, θ_(i) and φ_(i) are the angles of azimuth and elevation ofsource i (FIG. 1), A is the (N×P) matrix of the direction vectors of thesources a(θ_(i), φ_(i))(1≦i≦P), which contains in particular all of theinformation relating to the directions of arrival of the sources, A_(g)is the (N×P_(g)) sub-matrix of A relating to the gth group of sources,m_(g)(t) is the associated (P_(g)×1) subvector of m(t) and x_(g)(t)

A_(g) m_(g)(t). In particular, in the absence of coupling betweenantennas, the component n of the vector a(θ_(i), φ_(i)), denoteda_(n)(θ_(i), φ_(i)), can be written, in the case of an array whichutilizes only space diversity, in the following manner [5]

a _(n)(θ_(i),φ_(i))=exp{j2π[x _(n) cos(θ_(i))cos(φ_(i))+y _(n)sin(θ_(i))cos(φ_(i))+z _(n) sin(φ_(i))]/λ}  (3)

where λ is the wavelength and (x_(n), y_(n), z_(n)) are the coordinatesof antenna n.

The method according to the invention relies notably on the followingidea: the MUSIC-2q procedure (q≧1) implemented according to theinvention utilizes the information contained in the matrix of cumulantsof order 2q (N^(q)×N^(q)), C_(2q,x), whose coefficients are the circularcumulants of order 2q of the observations, Cum[x_(i) ₁ (t), . . . ,x_(i) _(q) (t), x_(i) _(q+1) (t)*, . . . , x_(i) _(2q)(t)*](1≦i_(j)≦N)(1≦j≦2q), where * corresponds to the complex conjugate.

These coefficients can be ranked in the matrix C_(2q,x) in various ways,which determines in particular the resolution and the processingcapacity of the MUSIC-2q procedure.

To parametrize these rankings, an arbitrary integer l is introduced suchthat (0≦l≦q) and the 2q-tuple, (i₁, . . . , i_(q), i_(q+1), . . . ,i_(2q)), with indices i_(j) (1≦j≦2q) is structured as 2 q-tuples indexedby l and defined respectively by (i₁, i₂, . . . , i_(l), i_(q+l), . . ., i_(2q−l)) and (i_(2q−l+1), . . . , i_(2q), i_(l+1), . . . , i_(q)).When the indices i_(j)(1≦j≦2q) vary from 1 to N, the previous 2 q-tuplestake N^(q) values. Numbering in a natural manner the N^(q) values ofeach of these 2 q-tuples by respectively the integers I_(l) and J_(l),such that 1≦I_(l), J_(l)≦N^(q), we obtain:

$\begin{matrix}{{I_{l} =};{\Delta \sum\limits_{j = 1}^{l}};;{{N^{q - j}\left( {i_{j} - 1} \right)} + \sum\limits_{j = 1}^{q - l}};;{{N^{q - l - j}\left( {i_{q + j} - 1} \right)} + 1}} & \left( {4a} \right) \\{{J_{l} =};{\Delta \sum\limits_{j = 1}^{l}};;{{N^{q - j}\left( {i_{{2q} - l + j} - 1} \right)} + \sum\limits_{j = 1}^{q - l}};;{{N^{q - l - j}\left( {i_{l + j} - 1} \right)} + 1}} & \left( {4b} \right)\end{matrix}$

Using the property of invariance under permutation of the cumulants, wededuce that Cum[x_(i) ₁ (t), . . . , x_(i) _(q) (t), x_(i) _(q+1) (t)*,. . . , x_(i) _(2q) (t)*]=Cum[x_(i) ₁ (t), . . . , x_(i) _(l) (t), x_(i)_(q+1) (t)*, . . . , x_(i) _(2q−l) (t)*, x_(i) _(2q−l+1) (t)*, . . . ,x_(i) _(2q) (t)*, x_(i) _(l+1) (t), . . . , x_(i) _(q) (t)] and assumingthat the latter quantity is the element [I_(l), J_(l)] of the matrixC_(2q,x), thus denoted C_(2q,x)(l), it is easy to show by using (2) thatthe matrix (N^(q)×N^(q)) C_(2q,x)(l) can be written in the followingmanner:

$\begin{matrix}{{{{C_{{2q},x}(l)} \approx \sum\limits_{g = 1}^{G}};};{{C_{{2q},x_{g}}(l)} + {\eta_{2}V\; {\delta \left( {q - 1} \right)}}}} & (5)\end{matrix}$

where η₂ is the average power of the noise per antenna, V is the (N×N)spatial coherence matrix of the noise such that Tr[V]=N, Tr[.] signifiesTrace, δ(.) is the Kronecker symbol and the (N^(q)×N^(q)) matrixC_(2q,x) _(g) (l) corresponds to the matrix of the circular cumulants oforder 2q of the vector x_(g)(t) for the ranking indexed by l, and may bewritten:

C _(2q,xg)(l)=[A _(g) ^({circle around (×)}l){circle around (×)}A_(g)*^({circle around (×)}(q−l)) ]C _(2q,m) _(g) (l)[A _(g)^({circle around (×)}l){circle around (×)}A_(g)*^({circle around (×)}(q−l))]^(†)  (6)

where C_(2q,m) _(g) (l) is the (P_(g) ^(q)×P_(g) ^(q)) matrix of thecircular cumulants of order 2q of m_(g)(t) for the ranking indexed by l,† corresponds to the conjugate transpose operation, {circle around (×)}is the Kronecker product and A_(g) ^({circle around (×)}l) is the(N^(l)×P_(g) ^(l)) matrix defined by A_(g) ^({circle around (×)}l)

A_(g){circle around (×)}A_(g){circle around (×)} . . . {circle around(×)}A_(g) with a number of Kronecker products equal to l−1.

In particular, for q=1 and l=1, the (N×N) matrix C_(2q,x)(l) correspondsto the covariance matrix of the observations (insofar as theobservations are centered) defined by

$\begin{matrix}{{R_{x} =};{{\Delta \; {C_{2,x}(1)}} = {{E\left\lbrack {{x(t)}{x(t)}^{\dagger}} \right\rbrack} \approx \sum\limits_{g = 1}^{G}}};;{{A_{g}{C_{2,m_{g}}(1)}A_{g}^{\dagger}} + {\eta_{2}V}}} & (7)\end{matrix}$

For q=2 and l=1, the (N²×N²) matrix C_(2q,x)(l) corresponds to theconventional quadricovariance matrix of the observations, defined by

$\begin{matrix}{{Q_{x} =};{{{\Delta \; {C_{4,x}(1)}} \approx \sum\limits_{g = 1}^{G}};;\left\lbrack {A_{g} \otimes A_{g}^{*}} \right\rbrack {{C_{4,m_{g}}(1)}\left\lbrack {A_{g} \otimes A_{g}^{*}} \right\rbrack}^{\dagger}}} & (8)\end{matrix}$

while for q=2 and l=2, the (N²×N²) matrix C_(2q,x)(l) corresponds to analternative expression for the quadricovariance matrix of theobservations, defined by

$\begin{matrix}{{{{\overset{\sim}{Q}}_{x};} =};{{\Delta \; {C_{4,x}(2)}} \approx \sum\limits_{g = 1}^{G}};;{\left\lbrack {A_{g} \otimes A_{g}} \right\rbrack {{C_{4,m_{g}}(2)}\left\lbrack {A_{g} \otimes A_{g}} \right\rbrack}^{\dagger}}} & (9)\end{matrix}$

Estimation

In practice, the statistics of order 2q of the observations, Cum[x_(i) ₁(t), . . . , x_(i) _(q) (t), x_(i) _(q+1) (t)*, . . . , x_(i) _(2q)(t)*], are unknown a priori and must be estimated on the basis of Lsamples of the observations, x(k)

x(kT_(e)), 1≦k≦L, where T_(e) is the sampling period.

For centered and stationary observations, using the ergodicity property,an asymptotically unbiased and consistent empirical estimator of thecumulants Cum[x_(i) ₁ (t), . . . , x_(i) _(q) (t), x_(i) _(q+1) (t)*, .. . , x_(i) _(2q) (t)*] can be constructed on the basis of the wellknown Leonov-Shiryaev formula [11], giving the expression for a cumulantof order n of x(t) as a function of its moments of order p (1≦p≦n), byreplacing all the moments in the latter formula by their empiricalestimate. More precisely, the Leonov-Shiryaev formula is given by:

$\begin{matrix}{{{{{{Cum}\left\lbrack {{x_{i_{1}}(t)}^{ɛ1},{x_{i_{2}}(t)}^{ɛ2},\ldots \;,{x_{i_{n}}(t)}^{ɛ\; n}} \right\rbrack} = \sum\limits_{p = 1}^{n}};};{\left( {- 1} \right)^{p - 1}{\left( {p - 1} \right)!}}}{{E\left\lbrack {{\prod\limits_{j \in {S\; 1}};};{x_{i_{j}}(t)}^{ɛ\; j}} \right\rbrack}{E\left\lbrack {{\prod\limits_{j \in {S\; 2}};};{x_{i_{j}}(t)}^{ɛ\; j}} \right\rbrack}\mspace{11mu} \ldots \mspace{11mu} {E\left\lbrack {\prod\limits_{j \in {Sp}}{x_{i_{j}}(t)}^{ɛ\; j}} \right\rbrack}}} & (10)\end{matrix}$

where (S1, S2, . . . , Sp) describe all the partitions into p sets of(1, 2, . . . , n), ε_(j)=±1(1≦p≦n) with the convention x¹=x and x⁻¹=x*and an empirical estimate of (10) is obtained by replacing all themoments E[x_(i) ₁ (t)^(ε1)x_(i) ₂ (t)^(δ2) . . . x_(i) _(p)(t)^(εp)](1≦p≦n) in (10) by their empirical estimate given by:

$\begin{matrix}{{{{{{\hat{E};}\left\lbrack {{x_{i_{1}}(t)}^{ɛ1}{x_{i_{2}}(t)}^{ɛ2}\mspace{11mu} \ldots \mspace{11mu} {x_{i_{p}}(t)}^{ɛ\; p}} \right\rbrack}(L)} = {{\text{;}{\Delta 1}\text{;}} - {\text{;}L\sum\limits_{k = 1}^{L}}}};};{{x_{i_{1}}(k)}^{ɛ1}{x_{i_{2}}(k)}^{ɛ2}\mspace{11mu} \ldots \mspace{11mu} {x_{i_{p}}(k)}^{ɛ\; p}}} & (11)\end{matrix}$

However, in radiocommunications, most sources are not stationary but arecyclostationary (digital modulations). For cyclostationary centeredobservations, the statistics matrix defined by (5) becomes timedependent, denoted C_(2q,x)(l)(t), and the method according to theinvention can also be implemented by considering that C_(2q,x)(l) is, inthis case, the time average, <C_(2q,x)(l)(t)>, over an interval ofinfinite horizon, of the instantaneous statistics matrix,C_(2q,x)(l)(t). Under these conditions, using the cyclo-ergodicityproperty, the matrix C_(2q,x)(l) must be estimated on the basis of theobservations by a nonempirical estimator such as that presented in [7]for q=2. This extension also applies to uncentered cyclostationarysources such as sources modulated with certain frequency modulations[9], provided that a nonempirical statistical estimator, such as thatpresented in [9] for q=1 and in [8] for q=2, is used.

The method according to the invention, or MUSIC-2q procedure comprisesvarious steps and variant embodiments detailed hereinafter.

The properties of the covariance matrix of order 2q, C_(2q,x)(l) arefirstly analyzed and the MUSIC-2q algorithm for the ranking indexed by lis deduced therefrom.

Assumptions

To develop the MUSIC-2q algorithm for the ranking indexed by l, severalassumptions presented hereinafter must be made:

H1: P_(g)<N, 1≦g≦G

H2: the matrix A_(g) ^({circle around (×)}l){circle around(×)}A_(g)*^({circle around (×)}(q−l)) is of full rank P_(g) ^(q), 1≦g≦G

${{H\; 3\text{:}\mspace{11mu} {P\left( {G,q} \right)}} =};{\Delta \sum\limits_{g = 1}^{G}};;{P_{g}^{q} < N^{q}}$

H4: the matrix [A₁ ^({circle around (×)}l){circle around(×)}A₁*^({circle around (×)}(q−l)), . . . , A_(G)^({circle around (×)}l){circle around(×)}A_(G)*^({circle around (×)}(q−l))] is of full rank P(G, q).

In particular, for statistically dependent sources (G=1), theassumptions H1 to H4 reduce to:

H1′: P<N

H2′: the matrix A^({circle around (×)}l){circle around(×)}A*^({circle around (×)}(q−l)) is of full rank P^(q)

while for statistically independent sources (G=P), the assumptions H1 toH4 reduce to:

H1″: P<N^(q)

H2″: the matrix [a₁ ^({circle around (×)}l){circle around(×)}a₁*^({circle around (×)}(q−l)), . . . , a_(P)^({circle around (×)}l){circle around(×)}a_(P)*^({circle around (×)}(q−l))] is of full rank P

Properties of C_(2q,x)(l)

The (P_(g) ^(q)×P_(g) ^(q)) matrix C_(2q,m) _(g) (l), which contains thecircular cumulants of order 2q of m_(g)(t) for the ranking indexed by l,is in general of full rank, P_(g) ^(q), insofar as the components ofm_(g)(t) are statistically dependent. Hence, using H1 and H2, the matrixC_(2q,xg)(l) for q>1 also has the rank P_(g) ^(q). Under theseconditions, we deduce from H4 that for q>1, the rank, r_(2q,x)(l), ofthe matrix C_(2q,x)(l) is equal to P(G, q) and such thatr_(2q,x)(l)<N^(q) according to H3. In particular, for sources which areall statistically dependent, r_(2q,x)(l)=P^(q) while for statisticallyindependent sources, r_(2q,x)(l)=P. Insofar as the matrix C_(2q,x)(l) isHermitian, but not positive definite, we deduce from the previousresults that the matrix C_(2q,x)(l) has P(G, q) nonzero eigenvalues andN^(q)−P(G, q) zero eigenvalues.

The MUSIC-2q Algorithm

To construct an algorithm of MUSIC type on the basis of the matrixC_(2q,x)(l), for q>1, we first calculate its eigenelement decomposition,given by

C _(2q,x)(l)=U _(2q,s)(l)Λ_(2q,s)(l)U _(2q,s)(l)^(†) +U_(2q,n)(l)Λ_(2q,n)(l)U _(2q,n)(l)^(†)  (12)

where Λ_(2q,s)(l) is the (P(G, q)×P(G, q)) diagonal matrix of thenonzero eigenvalues of C_(2q,x)(l), U_(2q,s)(l) is the (N^(q)×P(G, q))unit matrix of the eigenvectors of C_(2q,x)(l) associated with the P(G,q) nonzero eigenvalues of C_(2q,x)(l), Λ_(2q,n)(l) is the ((N^(q)−P(G,q))×(N^(q)−P(G, q))) diagonal matrix of the zero eigenvalues ofC_(2q,x)(l) and U_(2q,n)(l) is the (N^(q)×(N^(q)−P(G, q))) unit matrixof the eigenvectors of C_(2q,x)(l) associated with the (N^(q)−P(G, q))zero eigenvalues of C_(2q,x)(l). Insofar as C_(2q,x)(l) is Hermitian,all the columns of U_(2q,s)(l) are orthogonal to all the columns ofU_(2q,n)(l). Furthermore, insofar as Span {U_(2q,s)(l)}=Span {[A₁^({circle around (×)}l){circle around(×)}A₁*^({circle around (×)}(q−l)), . . . , A_(G)^({circle around (×)}l){circle around(×)}A_(G)*^({circle around (×)}(q−l))]}, all the columns of all thematrices A_(g) ^({circle around (×)}l){circle around(×)}A_(g)*^({circle around (×)}(q−l)), 1≦g≦G, are orthogonal to all thecolumns of U_(2q,n)(l). Let us denote by (θ_(ig), φ_(ig)) the directionof arrival of source i of group g. Under these conditions, the vectora(θ_(ig), φ_(ig))^({circle around (×)}l){circle around (×)}a(θ_(ig),φ_(ig))*^({circle around (×)}(q−l)) appears as column [(1−P_(g)^(q))(i−1)/(1−P_(g))+1] of A_(g) ^({circle around (×)}l){circle around(×)}A_(g)*^({circle around (×)}(q−l)). Hence, all the vectors {a(θ_(ig),φ_(ig))^({circle around (×)}l){circle around (×)}a(θ_(ig),φ_(ig))*^({circle around (×)}(q−l)), 1≦i≦P_(g), 1≦g≦G} are orthogonal tothe columns of U_(2q,n)(l) and are the only solutions of the followingequation:

[a(θ,φ)^({circle around (×)}l){circle around(×)}a(θ,φ)*^({circle around (×)}(q−l))]^(†) U _(2q,n)(l)U _(2q,n)(l)^(†)[a(θ,φ)^({circle around (×)}l){circle around(×)}a(θ,φ)*^({circle around (×)}(q−l))]=0  (13)

which corresponds to the core of the MUSIC-2q algorithm for the rankingindexed by l. In practice, the matrix U_(2q,s)(l) must be estimated onthe basis of the observations and the directions of arrival of thesources must be estimated by searching for the minima of the left-handside of equation (13). Insofar as, for space diversity arrays, the normof the vector [a(θ, φ)^({circle around (×)}l){circle around (×)}a(θ,φ)*^({circle around (×)}(q−l))] is independent of the direction ofarrival, the first term on the left-hand side of equation (13) can benormalized by the term [a(θ, φ)^({circle around (×)}l){circle around(×)}a(θ, φ)*^({circle around (×)}(q−l))]^(\ [a(θ, φ))^({circle around (×)}l){circle around (×)}a(θ,φ)*^({circle around (×)}(q−l))] without changing the results of thealgorithm.

The various steps of the MUSIC-2q algorithm for the ranking indexed by lare summarized hereinafter:

1. Estimation, C;̂_(2q,x)(l), of the matrix C_(2q,x)(l) on the basis of Lsample vectors x(k), 1≦k≦L, by using an appropriate estimator of thecircular cumulants of order 2q of the observations,

2. Decomposition of the matrix, C;̂_(2q,x)(l), into eigenelements, andextraction of an estimate, U;̂_(2q,n)(l), of the matrix U_(2q,n)(l). Thisstep may require an estimation of rank in the case where the number ofsources and/or their degree of statistical dependency is unknown apriori.

3. Calculation of the estimated pseudo-spectrum

$\begin{matrix}{{{\hat{P}\text{;}_{{Music}\text{-}2{q{(l)}}}\left( {\theta,\phi} \right)} =};{\Delta \frac{\begin{matrix}{\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q\; - l})}}}} \right\rbrack^{\dagger}{{\hat{U}}_{{2q},n}(l)}{{\hat{U}}_{{2q},n}(l)}^{\dagger}} \\\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}{\begin{matrix}\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}}} & (14)\end{matrix}$

for a given meshing of the space and search for the local minima. Thesearch for the minima included for example a local interpolation abouteach of the local minima.

In certain cases, the number of sources P is known, such that P<N, buttheir degree of statistical dependency is not known. Under theseconditions, P(G, q)≦P^(q) and a conservative approach for constructingU;̂_(2q,n)(l) is to use only the (N^(q)−P^(q)) eigenvectors ofC;̂_(2q,x)(l) that are associated with the smallest eigenvalues, therebyimplicitly returning to the assumption that all the sources arestatistically dependent.

The steps of the method according to the invention that were detailedpreviously for arbitrary non-Gaussian sources and space diversityarrays, can also be applied to the case of space, pattern and/orpolarization diversity arrays as well as to the case where the sourcespossess a rectilinearity property that is known a priori and utilized inthe algorithm.

Space, Pattern and/or Polarization Diversity Arrays

In the general case of an array comprising different antennas, that isto say of an array with pattern and/or polarization diversity inaddition to or instead of space diversity, the direction vector ofsource i depends not only on the direction of arrival, (θ_(i), φ_(i)),of source i but also its polarization, p_(i), where p_(i) is a pair ofparameters corresponding to the angles of polarization of source i inthe waveplane. Under these conditions, the direction vector of source i,denoted a(θ_(i), φ_(i), p_(i)), has a component n given by [5]

a_(n)(θ_(i),φ_(i),p_(i))=f_(n)(θ_(i),φ_(i),p_(i))exp{j2π[x _(n)cos(θ_(i))cos(φ_(i))+y_(n) sin(θ_(i))cos(φ_(i))+z_(n)sin(φ_(i))]/λ}  (15)

where f_(n)(θ_(i), φ_(i), p_(i)) is the complex response of antenna n toan electric field originating from the direction (θ_(i), φ_(i)) andhaving the polarization state p_(i).

Under these assumptions, the MUSIC-2q algorithm described in theprevious paragraph applies in the same manner to arrays utilizingpattern and polarization diversity, but with an estimatedpseudo-spectrum that depends on the polarization parameter p and isgiven by:

$\begin{matrix}{{{\hat{P}\text{;}_{{Music}\text{-}2{q{(l)}}}\left( {\theta,\phi,p} \right)} =};{\Delta \frac{\begin{matrix}{\left\lbrack {{a\left( {\theta,\phi,p} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi,p} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger}{{\hat{U}}_{{2q},n}(l)}{{\hat{U}}_{{2q},n}(l)}^{\dagger}} \\\left\lbrack {{a\left( {\theta,\phi,p} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi,p} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}{\begin{matrix}\left\lbrack {{a\left( {\theta,\phi,p} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi,p} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\\left\lbrack {{a\left( {\theta,\phi,p} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi,p} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}}} & (16)\end{matrix}$

In the case of an array utilizing pattern and polarization diversities,the pseudo-spectrum is normalized by the term

[a(θ,φ, p)^({circle around (×)}l){circle around (×)}a(θ,φ,p)*^({circle around (×)}(q−l))]^(†) [a(θ,φ, p)^({circle around (×)}l){circle around (×)}a(θ,φ,p)*^({circle around (×)}(q−l))]

so as to render the pseudo-spectrum constant with respect to directionsin the absence of sources.Hence, under the assumption of control of the responses of the antennasin terms of direction and polarization (for example by calibrating thearray in terms of direction and polarization), the search for the minimaof the pseudo-spectrum, P;̂_(2q-Music(l))(θ, φ, p), demands a meshing notonly of the direction space (θ, φ) but also of the polarization space p,followed by an interpolation in the neighborhoods of the minima, therebyrepresenting a considerable cost. One way to circumvent the meshing ofthe polarization space is to utilize the fact that a polarization wave pcan be decomposed into the weighted sum of two waves with orthogonalpolarizations, p₁ and P₂. Under these conditions, the direction vectora(θ, φ, p) becomes a linear combination of the direction vectors a(θ, φ,p₁) and a(θ, φ, p₂) relating to these two polarizations such that:

a(θ,φ,p)=α₁ a(θ,φ,p ₁)+α₂ a(θ,φ,p ₂)=A ₁₂(θ,φ)α(p)  (17)

where α₁ and α₂, such that α₁ ²+α₂ ²=1, are the real coefficients of thedecomposition of the incident field into a sum of two orthogonallypolarized fields, α(p) is the (2×1) vector of components α₁ and α₂ andA₁₂(θ, φ) is the (N×2) matrix defined by A₁₂(θ, φ)

[a(θ, φ, p₁), a(θ, φ, p₂)].

Using the fact that (XY){circle around (×)}(ZW)=(X{circle around(×)}Z)(Y{circle around (×)}W), we deduce from (17) that

$\begin{matrix}{{{{a\left( {\theta,\phi,p} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi,p} \right)}^{\otimes {({q - l})}}} = \left\lbrack {{A_{12}\left( {\theta,\phi} \right)}^{\otimes l} \otimes {A_{12}\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack}{\quad\mspace{79mu} {{\left\lbrack {{\alpha (p)}^{\otimes l} \otimes {\alpha (p)}^{*{\otimes {({q - l})}}}} \right\rbrack =};{\Delta \; {A_{12,q,l}\left( {\theta,\phi} \right)}{\alpha_{q,l}(p)}}}}} & (18)\end{matrix}$

where A_(12,q,l)(θ, φ) is the (N^(q)×2^(q)) matrix defined byA_(12,q,l)(θ, φ)

[A₁₂(θ, φ)^({circle around (×)}l){circle around (×)}A₁₂(θ,φ)*^({circle around (×)}(q−l))] and where α_(q,l)(p) is the (2^(q)×1)vector defined by α_(q,l)(p)

[α(p)^({circle around (×)}l){circle around(×)}α(p)*^({circle around (×)}(q−l))]. Under these conditions, thepseudo-spectrum (16) takes the form:

$\begin{matrix}\begin{matrix}{{{\hat{P}\text{;}_{{Music}\text{-}2{q{(l)}}}\left( {\theta,\phi,p} \right)} =};\Delta} & \; \\\alpha_{q} & \; \\{{l(p)}^{\dagger}A_{12}} & \alpha_{q} \\q & {{l(p)}^{\dagger}A_{12}} \\{{l\left( {\theta,\phi} \right)}^{\dagger}\hat{U}\text{;}_{2q}} & q \\{{{n(l)}\hat{U}\text{;}_{2q}};} & {\text{;}_{l}\left( {\theta,\phi} \right)^{\dagger}A_{12}} \\{{n(l)}^{\dagger}A_{12}} & q \\q & {{l\left( {\theta,\phi} \right)}\alpha_{q}} \\{{l\left( {\theta,\phi} \right)}\alpha_{q}} & {l(p)} \\{l(p)} & \;\end{matrix} & (19)\end{matrix}$

and the use in (19) of the vector α_(q,l)(p) which minimizes (19)culminates in a pseudo-spectrum that now depends only on the directionof arrival and is given by:

P;̂ _(2q-Music(l))(θ,φ)=λ_(min)[(A _(12,q,l)(θ,φ)^(†) A _(12,q,l)(θ,φ))⁻¹A _(12,q,l)(θ,φ)^(†) U;̂ _(2q,n)(l)U;̂ _(2q,n)(l)^(†) A_(12,q,l)(θ,φ)]  (20)

where λ_(min)[X] corresponds to the minimum eigenvalue of the matrix X.The estimated directions of the sources, (θ;̂_(i), φ;̂_(i))(1≦i≦P),correspond to the local minima of (20), which have to be searched for bycontrolling the direction vectors a(θ, φ, p₁) and a(θ, φ, p₂) in the twopolarizations p₁ and p₂ (by calibration for example), meshing thedirection of arrival space followed by an interpolation about theminima.

An estimation, α;̂_(q,l)(p_(i)), of the vector α_(q,l)(p_(i)) associatedwith source i is then given by the eigenvector of the matrix(A_(12,q,l)(θ;̂_(i), φ;̂_(i))^(†)A_(12,q,l)(θ;̂_(i),φ;̂_(i)))⁻¹A_(12,q,l)(θ;̂_(i), φ;̂_(i))^(†)U;̂_(2q,n)(l)U;̂_(2q,n)(l)^(†)A_(12,q,l)(θ;̂_(i), φ;̂_(i)) that is associated with theminimum eigenvalue. By decomposing the vector α;̂_(q,l)(p_(i)) into2^(q−2)(4×1) subvectors, α;̂_(q,l,s)(p_(i)), 1≦s≦q−2, we obtainα;̂_(q,l)(p_(i))=[α;̂_(q,l,1)(p_(i))^(T), . . . ,α;̂_(q,l,(q−2))(p_(i))^(T)]^(T). It is then easy to verify that thesubvector α;̂_(q,l,s)(p_(i)) is an estimate of the vectorα_(q,l,s)(p_(i)), which is proportional to:

α(p _(i)){circle around (×)}α(p _(i)) if q−l=0  (21a)

α(p _(i)){circle around (×)}α(p _(i))* if q−l=1  (21b)

α(p _(i))*{circle around (×)}α(p _(i))* if q−l>1  (21c)

Under these conditions, by arranging the components of the vectorα;̂_(q,l,s)(p_(i)) in a (2×2) matrix Γ;̂_(q,l,s)(p_(i)) such thatΓ;̂_(q,l,s)(p_(i))[k, j]=α;̂_(q,l,s)(p_(i))[2(k−1)+j], whereΓ;̂_(q,l,s)(p_(i))[k, j] and α;̂_(q,l,s)(p_(i))[k] are the element [k, j]and the component k respectively of the matrix Γ;̂_(q,l,s)(p_(i)) and ofthe vector α;̂_(q,l,s)(p_(i)), the matrix Γ;̂_(q,l,s)(p_(i)) becomes, towithin a scalar, an estimate of the matrix Γ_(q,l,s)(p_(i)) defined by:

Γ_(q,l,s)(p _(i))=α(p _(i))α(p _(i))^(T) if q−l=0  (22a)

Γ_(q,l,s)(p _(i))=α(p _(i))α(p _(i))^(†) if q−l=1  (22b)

Γ_(q,l,s)(p _(i))=α(p _(i))*α(p _(i))^(†) if q−l>1  (22c)

It then suffices to jointly diagonalize the matrices Δ;̂_(q,l,s)(p_(i)),1≦s≦q−2, to obtain an estimate, α;̂(p_(i)), of α(p_(i)) which correspondsto the eigenvector associated with the maximum eigenvalue, whereΔ;̂_(q,l,s)(p_(i)) is defined by:

Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))Γ;̂_(q,l,s)(p _(i))^(†) ifq−l=0  (23a)

Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))Γ;̂_(q,l,s)(p _(i))^(†) ifq−l=1  (23b)

Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))*Γ;̂_(q,l,s)(p _(i))^(T) ifq−l>1  (23c)

Case of Rectilinear Sources

For noncircular observations to order 2q, we utilize the(2^(q)N^(q)×2^(q)N^(q)) matrix, C_(2q,xe), of cumulants of order 2q ofthe extended observation vector (2N×1) x_(e)(t)

[x(t)^(T), x(t)^(†)]^(T), whose coefficients are the circular cumulantsof order 2q of the vector x_(e)(t), Cum[x_(ei) ₁ (t), . . . , x_(ei)_(q) (t), x_(ei) _(q+1) (t)*, . . . , x_(ei) _(2q)(t)*](1≦i_(j)≦2N)(1≦j≦2q). These cumulants of order 2q can be ranked invarious ways, indexed by the integer l, in the matrix C_(2q,xe), whichis denoted, for the ranking indexed by l, C_(2q,xe)(l).

A particular family of noncircular sources to order 2q corresponds torectilinear sources, characterized by the fact that their complexenvelope is real, this being in particular the case for amplitudemodulated (AM, ASK) sources or 2-state phase modulated (BPSK) sources.For P sources all rectilinear, under the assumptions stated previously,the expression for the vector x_(e)(t) is deduced from (2) and is givenby:

$\begin{matrix}\begin{matrix}{{{{x_{e}(t)} \approx \sum\limits_{i = 1}^{P}};};{{{{m_{i}(t)}{a_{e}\left( {\theta_{i},\phi_{i},\varphi_{i}} \right)}} + {v_{e}(t)}} =};{{\Delta \; A_{e}{m(t)}} + {v_{e}(t)}}} \\{{{= \sum\limits_{g = 1}^{G}};};{{{A_{eg}{m_{g}(t)}} + {v_{e}(t)}} = \sum\limits_{g = 1}^{G}};;{{x_{eg}(t)} + {v_{e}(t)}}}\end{matrix} & (24)\end{matrix}$

where ν_(e)(t)

[ν(t)^(T), ν(t)^(†)]^(T), a_(e)(θ_(i), φ_(i), φ_(i))

[e^(jφ) _(i) a(θ_(i), φ_(i))^(T), e^(−jφ) _(i) a(θ_(i), φ_(i))^(†)]^(T),φ_(i) is the phase of the propagation channel associated with source i,A_(e) is the (2N×P) matrix of the vectors a_(e)(θ_(i), φ_(i),φ_(i))(1≦i≦P), A_(eg) is the (2N×P_(g)) sub-matrix of A_(e) relating tothe gth group of sources and x_(eg)(t)

A_(eg)m_(g)(t). The model (24) is similar to the model (2), in which theextended (2N×1) direction vectors, a_(e)(θ_(i), φ_(i), φ_(i)), havereplaced the (N×1) direction vectors, a(θ_(i), φ_(i)). In this context,under conditions C1 to C4 presented hereinafter:

C1: P_(g)<2N, 1≦g≦G

C2: the matrix A_(eg) ^({circle around (×)}l){circle around(×)}A_(eg)*^({circle around (×)}(q−l)) is of full rank P_(g) ^(q), 1≦g≦G

${{C\; 3\text{:}\mspace{11mu} {P\left( {G,q} \right)}} =};{\Delta \sum\limits_{g = 1}^{G}};;{P_{g}^{q} < {2^{q}N^{q}}}$

C4: the matrix [A_(e1) ^({circle around (×)}l){circle around(×)}A_(e1)*^({circle around (×)}(q−l)), . . . , A_(eG)^({circle around (×)}l){circle around(×)}A_(eG)*^({circle around (×)}(q−l))] is of full rank P(G, q) and byusing similar reasoning to that set out above regarding the propertiesof the matrix and the Music-2q algorithm, but on the basis of the matrixC_(2q,xe)(l), we culminate in the MUSIC-2q algorithm for rectilinearsources, denoted MUSIC-2q-REC, and for the ranking indexed by l, whosesteps are summarized hereinafter:

1. Estimation, C;̂_(2q,xe)(l), of the matrix C_(2q,xe)(l) on the basis ofL sample vectors x(k), 1≦k≦L, by using an appropriate estimator of thecircular cumulants of order 2q of the observations.

2. Decomposition of the matrix, C;̂_(2q,xe)(l), into eigenelements, andextraction of an estimate, U;̂_(2q,en)(l), of the matrix U_(2q,en)(l) ofthe eigenvectors of the noise space of the matrix C_(2q,xe)(l). Thisstep may require an estimation of rank in the case where the number ofsources and/or their degree of statistical dependency is unknown apriori.

3. Calculation of the estimated pseudo-spectrum

$\begin{matrix}{{{\hat{P}}_{{Music}\text{-}2q\text{-}{{rec}{(l)}}}\left( {\theta,\phi,\varphi} \right)}\overset{\Delta}{=}\frac{\begin{matrix}\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\{{{\hat{U}}_{{2q},{en}}(l)}{{\hat{U}}_{{2q},{en}}(l)}^{\dagger}} \\\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}{\begin{matrix}\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}} & (25)\end{matrix}$

for a given meshing of the direction space (θ, φ) and of the phase spaceφ and search for the local minima (including a local interpolation abouteach of the local minima).

According to a variant embodiment, one way to circumvent the meshing ofthe phase space is to utilize the fact that the extended directionvector α_(ε)(θ, φ, φ) can be written:

a _(e)(θ,φ,φ)=A;{tilde over ( )} _(e)(θ,φ)β(φ)  (26)

where β(φ) is the (2×1) vector defined by β(φ)

[e^(jφ), e^(−jφ)]^(T) and A;{tilde over ( )}_(e)(θ, φ) is the (2N×2)matrix defined by

A;{tilde over ( )} _(e)(θ,φ)

(;a(θ,φ)0;;0a(θ,φ)*;)  (27)

where 0 is the zero vector of dimension (N×1). Using the fact that(XY){circle around (×)}(ZW)=(X{circle around (×)}Z) (Y{circle around(×)}W), we deduce from (26) that:

a _(e)(θ,φ,φ)^({circle around (×)}l) {circle around (×)}a_(e)(θ,φ,φ)*^({circle around (×)}(q−l)) =[A;{tilde over ()}_(e)(θ,φ)^({circle around (×)}l) {circle around (×)}A;{tilde over ()}_(e)(θ,φ)*^({circle around (×)}(q−l))][β(φ)^({circle around (×)}l){circlearound (×)}β(φ)*^({circle around (×)}(q−l)) ]

A;{tilde over ( )}_(e,q,l)(θ,φ)β_(q,l)(φ)  (28)

where A;{tilde over ( )}_(e,q,l)(θ, φ) is the (2^(q)N^(q)×2^(q)) matrixdefined by A;{tilde over ( )}_(e,q,l)(θ, φ)

[A;{tilde over ( )}_(e)(θ, φ)^({circle around (×)}l){circle around(×)}A;{tilde over ( )}_(e)(θ, φ)*^({circle around (×)}(q−l))] and whereβ_(q,l)(φ) is the (2^(q)×1) vector defined by β_(q,l)(φ)

[β(φ)^({circle around (×)}l){circle around(×)}β(φ)*^({circle around (×)}(q−l))]. Under these conditions, by usingthe reasoning of the previous paragraph on phase diversity arrays, thepseudo-spectrum (25) takes the form:

P;̂ _(Music-2q-rec(l))(θ,φ)=λ_(min)[(A;{tilde over ( )} _(e,q,l)(θ,φ)^(†)A;{tilde over ( )} _(e,q,l)(θ,φ))⁻¹ A;{tilde over ( )} _(e,q,l)(θ,φ)^(†)U;̂ _(2q,en)(l)U;̂ _(2q,en)(l)^(†) A;{tilde over ( )} _(e,q,l)(θ,φ)]  (29)

The estimated directions of the sources, (θ;̂_(i), φ;̂_(i))(1≦i≦P),correspond to the local minima of (29), which have to be searched for bycontrolling the direction vector a(θ, φ) and meshing the direction ofarrival space solely, followed by an interpolation about the minima.

An estimation, β;̂_(q,l)(φ_(i)), of the vector β_(q,l)(φ_(i)) associatedwith source i is then given by the eigenvector of the matrix (A;{tildeover ( )}_(e,q,l)(θ;̂_(i), φ;̂_(i))^(†)A;{tilde over ( )}_(e,q,l)(θ;̂_(i),φ;̂_(i)))⁻¹A;{tilde over ( )}_(e,q,l)(θ;̂_(i), φ;̂_(i))^(†)U;̂_(2q,en)(l)U;̂_(2q,en)(l)^(†)A;{tilde over ( )}_(e,q,l)(θ;̂_(i), φ;̂_(i)) that isassociated with the minimum eigenvalue.

The MUSIC-2q-REC algorithm can also be implemented on the basis ofpolarization diversity arrays on the basis of meshing just the directionspace.

Identifiability

It can be shown that the order 2q antenna processing problem for Pnon-Gaussian and statistically independent sources on the basis of anarray of N antennas with coordinates (x_(n), y_(n), z_(n)) and responsesf_(n)(θ, φ, p), 1≦n≦N, is, for the ranking indexed by l, C_(2q,x)(l),similar to a second-order antenna processing problem for which these Pstatistically independent sources are incident, with a virtual powerc_(2q,m) _(i) (1≦i≦P), on a virtual array of N^(q) VAs with coordinates(x_(k) ₁ _(k) _(2 . . .) _(k) _(q) ^(l), y_(k) ₁ _(k) _(2 . . .) _(k)_(q) ^(l), z_(k) ₁ _(k) _(2 . . .) _(k) _(q) ^(l)) and responses f_(k) ₁_(k) _(2 . . .) _(k) _(q) ^(l)(θ, φ, p), 1≦k_(j)≦N for 1≦j≦q, of whichin general N_(2q);^(l) are different, defined respectively by:

$\begin{matrix}{\left( {x_{k_{1}k_{2}\mspace{11mu} \ldots \mspace{11mu} k_{q}}^{l},y_{k_{1}k_{2}\mspace{11mu} \ldots \mspace{11mu} k_{q}}^{l},z_{k_{1}k_{2}\mspace{11mu} \ldots \mspace{11mu} k_{q}}^{l}} \right) = \left( {{{\sum\limits_{j = 1}^{l};};{x_{k_{j}} - \sum\limits_{u = 1}^{q - l}};;x_{k_{l + u}}},{{\sum\limits_{j = 1}^{l};};{y_{k_{j}} - \sum\limits_{u = 1}^{q - l}};;y_{k_{l + u}}},{{\sum\limits_{j = 1}^{l};};{z_{k_{j}} - \sum\limits_{u = 1}^{l}};;z_{k_{l + u}}}} \right)} & (30) \\{{{{f_{k_{1}k_{2}\mspace{11mu} \ldots \mspace{11mu} k_{q}}^{l}\left( {\theta,\phi,p} \right)} = \prod\limits_{j = 1}^{l}};};{\prod\limits_{u = 1}^{q - l};};{{f_{k_{j}}\left( {\theta,\phi,p} \right)}{f_{k_{l + u}}\left( {\theta,\phi,p} \right)}^{*}}} & (31)\end{matrix}$

Thus antenna processing to higher orders can be used to replace antennasand hardware and thus decrease the overall cost of the system.

To illustrate this notion, FIG. 2 presents the virtual array associatedwith a uniform circular array of 5 antennas of radius R=0.8λ, jointlywith the order of multiplicity of the virtual antennas (that is to sayof the number of virtual sensors corresponding to this antenna) for (q,l)=(3, 2). In this typical case N_(2q);^(l)=55.

Utilization of C_(2q,xe)(l) and Rectilinear Sources

The previously mentioned theory of virtual arrays applies to the problemof the antenna processing to order 2q for rectilinear sources utilizingthe information contained in the matrix, C_(2q,xe)(l), of all thecumulants of order 2q of the observations for the ranking indexed by lby replacing the initial N-antenna array with the 2N-antenna virtualarray associated with the extended observations. We then denote byN_(2q);^(l) _(,e) the number of different VAs of the virtual array oforder 2q for the ranking indexed by l associated with the antennaprocessing problem for rectilinear sources utilizing all the cumulantsof order 2q of the observations.

Performance of the MUSIC-2q Algorithm

We deduce from the previous results that the number of non-Gaussian andstatistically independent sources that can be processed by the MUSIC-2qalgorithm for the ranking indexed by l is Inf(N_(2q);^(l), N^(q)−1).

Tables 1 and 2 present an upper bound, N_(max)[2q, l], of N_(2q);^(l),for polarization diversity and space diversity arrays respectively,attained for arrays not having any particular symmetry.

Performance of the MUSIC-2q-REC Algorithm

Similar reasoning to the above shows that the number of non-Gaussian andstatistically independent sources that can be processed by theMUSIC-2q-REC algorithm for the ranking indexed by l is Inf(N_(2q);^(l)_(,e), 2^(q)N^(q)−1).

TABLE 1 N_(max)[2q, l] as a function of N for several values of q and land for space and polarization diversity arrays. m = 2q l N/_(max)[2q,l] 4 2 N(N + 1)/2 (q = 2) 1 N² 6 3 N!/[6(N − 3)!] + N(N − 1) + N (q = 3)2 N!/[2(N − 3)!] + 2N(N − 1) + N 8 4 N!/[24(N − 4)!] + N!/[2(N − 3)!] +1.5N(N − 1) + N (q = 4) 3 N!/[6(N − 4)!] + 1.5N!/(N − 3)!] + 3N(N − 1) +N 2 N!/[4(N − 4)!] + 2N!/(N − 3)!] + 3.5N(N − 1) + N

TABLE 1 N_(max)[2q, l] as a function of N for several values of q and land for space diversity arrays. m = 2q l N/_(max)[2q, l] 4 2 N(N + 1)/2(q = 2) 1 N² − N + 1 6 3 N!/[6(N − 3)!] + N(N − 1) + N (q = 3) 2 N!/[2(N− 3)!] + N(N − 1) + N 8 4 N!/[24(N − 4)!] + N!/[2(N − 3)!] + 1.5N(N− 1) + N (q = 4) 3 N!/[6(N − 4)!] + N!/(N − 3)! + 1.5N(N − 1) + N 2N!/[4(N − 4)!] + N!/(N − 3)! + 2N(N − 1) + 1

The invention relates notably to a device for high-resolution directionfinding to an arbitrary even order, 2q, with q>2 disposed in an arraycomprising N narrowband antennas each receiving the contribution from Psources, the device comprising a processor suitable for executing thesteps of the method described above.

It will be readily seen by one of ordinary skill in the art that thepresent invention fulfils all of the objects set forth above. Afterreading the foregoing specification, one of ordinary skill in the artwill be able to affect various changes, substitutions of equivalents andvarious aspects of the invention as broadly disclosed herein. It istherefore intended that the protection granted hereon be limited only bydefinition contained in the appended claims and equivalents thereof.

REFERENCES

-   [1] G. BIENVENU, L. KOPP, “optimality of high resolution array    processing using the eigensystem approach”, IEEE Trans. Acou. Speech    and Sign. Proc., Vol. 31, No. 5, pp. 1235-1247, October 1983.-   [2] J. P. BURG, “The relationship between maximum entropy spectra    and maximum likelihood spectra”, Geophysics, Vol. 37, No. 2, pp.    375-376, April 1972.-   [3] J. CAPON, “High resolution frequency-wavenumber spectrum    analysis”, Proc. IEEE, Vol. 57, No. 8, pp. 1408-1418, February 1969.-   [4] P. CHEVALIER, A. FERREOL, “On the virtual array concept for the    fourth-order direction finding problem”, IEEE Trans. Signal    Processing, Vol. 47, No. 9, pp. 2592-2595, September 1999-   [5] R. T. COMPTON, JR., “Adaptive Antennas—Concepts and    Performance”, Prentice Hall, Englewood Cliffs, N.J., 07632, 1988.-   [6] M. C. DOGAN, J. M. MENDEL, “Applications of cumulants to array    processing—Part I: Aperture extension and array calibration”, IEEE    Trans. Signal Processing, Vol. 43, No. 5, pp. 1200-1216, May 1995.-   [7] A. FERREOL, P. CHEVALIER, “On the behavior of current second and    higher order blind source separation methods for cyclostationary    sources”, IEEE Trans. Signal Processing, Vol. 48, No. 6, pp.    1712-1725, June 2000. Errata Vol. 50, No. 4, p 990, April 2002.-   [8] A. FERREOL, P. CHEVALIER, L. ALBERA, “Higher order blind    separation of non zero-mean cyclostationary sources”, Proc. EUSIPCO    02, Toulouse, (France), pp. 103-106, September 2002.-   [9] A. FERREOL, P. CHEVALIER, L. ALBERA, “Second order blind    separation of first and second order cyclostationary    sources—Application to AM, FSK, CPFSK and Deterministic sources”,    IEEE Trans. Signal Processing, Vol. 52, No. 4, pp. 845-861, April    2004.-   [10] B. FRIEDLANDER, “A sensitivity analysis of the MUSIC    algorithm”, IEEE Trans. Acou. Speech. and Signal Processing, Vol.    38, No. 10, pp. 1740-1751, October 1990.-   [11] P. McCULLAGH, “Tensor methods in Statistics”, Chapman and Hall,    Monographs on Statistics and applied Probability, 1987.-   [12] A. PAULRAJ, T. KAILATH, “Eigenstructure methods for direction    of arrival estimation in the presence of unknown noise field”, IEEE    Trans. Acou. Speech and Sign. Proc., Vol. 34, No. 1, pp. 13-20,    February 1986.-   [13] B. PORAT, B. FRIEDLANDER, “Direction finding algorithms based    on higher order statistics”, IEEE Trans. Signal Processing, Vol. 39,    No. 9, pp. 2016-2024, September 1991.-   [14] R. O. SCHMIDT, “Multiple emitter location and signal parameter    estimation”, IEEE Trans. Ant. Prop., Vol. 34, No. 3, pp. 276-280,    March 1986.-   [15] A. L. SWINDLEHURST, T. KAILATH, “A performance analysis of    subspaced-based methods in the presence of model errors, Part I: The    MUSIC algorithm”, IEEE Trans. Signal Processing, Vol. 40, No. 3, pp.    1758-1773, July 1992.

1. A method of high-resolution direction finding to an arbitrary evenorder, 2q (q>2), for an array comprising N narrowband antennas eachreceiving the contribution from P sources, wherein the algebraicproperties of a matrix of cumulants of order 2q, C_(2q,x)(l), whosecoefficients are the circular cumulants of order 2q, Cum[x_(i) ₁ (t), .. . , x_(i) _(q) (t), x_(i) _(q+1) (t)*, . . . , x_(i) _(2q) (t)*], ofthe observations received on each antenna, for cumulant rankings indexedby l, are utilized to define a signal subspace and a noise subspace. 2.The method as claimed in claim 1, comprising the following steps:estimating C;̂_(2q,x)(l), the cumulant matrix of order 2q, C_(2q,x)(l),on the basis of L sample vectors x(k), 1≦k≦L, by using an estimator ofthe circular cumulants of order 2q of the observations, decomposing theestimated matrix C;̂_(2q,x)(l) into eigenelements, and extracting anestimate, U;̂_(2q,n)(l), of the unit matrix U_(2q,n)(l) of theeigenvectors of C_(2q,x)(l), calculating the estimated pseudo-spectrumfor a given meshing of the direction space (θ, φ).${{\hat{P}}_{{Music}\text{-}2{q{(l)}}}\left( {\theta,\phi} \right)}\overset{\Delta}{=}\frac{\begin{matrix}{\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger}{{\hat{U}}_{{2q},n}(l)}{{\hat{U}}_{{2q},n}(l)}^{\dagger}} \\\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}{\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{\otimes {({q - l})}}} \right\rbrack^{\dagger}\left\lbrack {{a\left( {\theta,\phi} \right)}^{\otimes l} \otimes {a\left( {\theta,\phi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack}$3. The method as claimed in claim 2, wherein for a polarizationdiversity array, the estimated pseudo-spectrum is calculated for a givenmeshing of the direction space (θ, φ) and of the polarization space pand local minima are searched for.
 4. The method as claimed in claim 1,wherein for polarization diversity arrays, the method comprises at leastthe following steps: estimating C;̂_(2q,x)(l), the covariance matrixC_(2q,x)(l) on the basis of L sample vectors x(k), 1≦k≦L, by using anestimator of the circular cumulants of order 2q of the observations,decomposing the estimated matrix C;̂_(2q,x)(l) into eigenelements, andextracting an estimate, U;̂_(2q,n)(l), of the unit matrix U_(2q,n)(l) ofthe eigenvectors of C_(2q,x)(l), calculating the estimatedpseudo-spectrumP;̂ _(2q-Music(l))(θ,φ)=λ_(min)[(A _(12,q,l)(θ,φ)^(†) A _(12,q,l)(θ,φ))⁻¹A _(12,q,l)(θ,φ)^(†) U;̂ _(2q,n)(l)U;̂ _(2q,n)(l)^(†) A _(12,q,l)(θ,φ)]for a given meshing of the direction space (θ, φ), by controlling thedirection vectors a(θ, φ, p₁) and a(θ,φ, p₂) in the two polarizations p₁and p₂ and searching for local minima, where A₁₂(θ, φ) is the (N×2)matrix defined by A₁₂(θ, φ)

[a(θ, φ, p₁), a(θ, φ, p₂)] and where λ_(min)[X] corresponds to theminimum eigenvalue of the matrix X.
 5. The method as claimed in claim 4,comprising the following steps: estimating α;̂_(q,l)(p_(i)), the vectorα_(q,l)(p_(i)) associated with source i is then given by the eigenvectorof the matrix (A_(12,q,l)(θ;̂_(i), φ;̂_(i))^(†) A_(12,q,l)(θ;̂_(i),φ;̂_(i)))⁻¹ A_(12,q,l)(θ;̂_(i), φ;̂_(i))^(†) U;̂_(2q,n)(l)U;̂_(2q,n)(l)^(†)A_(12,q,l)(θ;̂_(i), φ;̂_(i)) that is associated with theminimum eigenvalue, decomposing the vector α;̂_(q,l)(p_(i)) into2^(q−2)(4×1) subvectors, α;̂_(q,l,s)(p_(i)), 1≦s≦q−2, such thatα;̂_(q,l)(p_(i))=[α;̂_(q,l,1)(p_(i))^(T), . . . ,α;̂_(q,l,(q−2))(p_(i))^(T)]^(T) and ranking the components of each vectorα;̂_(q,l,s)(p_(i)) in a (2×2) matrix Γ;̂_(q,l,s)(p_(i)) such that:Γ;̂_(q,l,s)(p_(i))[k, j]=α;̂_(q,l,s)(p_(i))[2(k−1)+j], whereΓ;̂_(q,l,s)(p_(i))[k, j] and α;̂_(q,l,s)(p_(i))[k] are the element [k, j]and the component k respectively of the matrix Γ;̂_(q,l,s)(p_(i)) and ofthe vector α;̂_(q,l,s)(p_(i)). α(p) being made up of the realcoefficients of the decomposition of the incident field into a sum oftwo orthogonally polarized fields, jointly diagonalizing the matricesΔ;̂_(q,l,s)(p_(i)), 1≦s≦q−2, to obtain an estimate, α;̂(p_(i)), of thevector α(p_(i)) of the polarization parameters of source i, estimatecorresponding to the eigenvector associated with the maximum eigenvalue,where Δ;̂_(q,l,s)(p_(i)) is defined by:Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))Γ;̂_(q,l,s)(p _(i))^(†) if q−l=0Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))Γ;̂_(q,l,s)(p _(i))^(†) if q−l=1Δ;̂_(q,l,s)(p _(i))=Γ;̂_(q,l,s)(p _(i))*Γ;̂_(q,l,s)(p _(i))^(T) if q−l>1 6.The method as claimed in claim 1 for rectilinear sources, wherein itcomprises at least the following steps: estimating C;̂_(2q,xe)(l), thematrix C_(2q,xe)(l) of cumulants of order 2q of the extendedobservations, x_(e)(t)=;^(Δ)[x(t)^(T), x(t)^(†)]^(T), where x(t) is theobservation vector on the basis of L sample vectors x(k), 1≦k≦L, byusing an estimator of the cumulants of order 2q of the observations,decomposing the matrix, C;̂_(2q,xe)(l), into eigenelements, andextracting an estimate, U;̂_(2q,en)(l), of the matrix U_(2q,en)(l) of theeigenvectors of the noise space of the matrix C_(2q,xe)(l), calculatingthe estimated pseudo-spectrum${{\hat{P}}_{{Music}\text{-}2q\text{-}{{rec}{(l)}}}\left( {\theta,\phi,\varphi} \right)}\overset{\Delta}{=}\frac{\begin{matrix}\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\{{{\hat{U}}_{{2q},{en}}(l)}{{\hat{U}}_{{2q},{en}}(l)}^{\dagger}} \\\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack\end{matrix}}{\begin{matrix}\left\lbrack {{a_{e}\left( {\theta,\phi,\varphi} \right)}^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}} \right\rbrack^{\dagger} \\\left\lbrack {a_{e}{\left( {\theta,\phi,\varphi} \right)^{\otimes l} \otimes {a_{e}\left( {\theta,\phi,\varphi} \right)}^{*{\otimes {({q - l})}}}}} \right\rbrack\end{matrix}}$ for a given meshing of the direction space (θ, φ) and ofthe phase space φ and searching for local minima.
 7. The method asclaimed in claim 1, comprising the following steps: estimating,C;̂_(2q,xe)(l), the matrix C_(2q,xe)(l) of the set of the cumulants oforder 2q of the observations, on the basis of L sample vectors x(k),1≦k≦L, by using an appropriate estimator of the cumulants of order 2q ofthe observations, decomposing the matrix, C;̂_(2q,xe)(l), intoeigenelements, and extracting an estimate, U;̂_(2q,en)(l), of the matrixU_(2q,en)(l) of the eigenvectors of the noise space of the matrixC_(2q,xe)(l), calculating the estimated pseudo-spectrumP;̂ _(Music-2q-rec(l))(θ,φ)=λ_(min)[(A;{tilde over ( )} _(e,q,l)(θ,φ)^(†)A;{tilde over ( )} _(e,q,l)(θ,φ))⁻¹ A;{tilde over ( )} _(e,q,l)(θ,φ)^(†)U;̂ _(2q,en)(l)U;̂ _(2q,en)(l)^(†) A;{tilde over ( )} _(e,q,l)(θ,φ)] whereA;{tilde over ( )}_(e,q,l)(θ, φ) is the matrix (2^(q)N^(q)×2^(q))defined byA;{tilde over ( )} _(e,q,l)(θ,φ)

[A;{tilde over ( )} _(e)(θ,φ)^({circle around (×)}l){circle around(×)}A;{tilde over ( )} _(e)(θ,φ)*^({circle around (×)}(q−l)]) whereA;{tilde over ( )}_(e)(θ, φ) is the matrix (2N×2) defined byA;{tilde over ( )} _(e)(θ,φ)

(;a(θ,φ)0;;0a(θ,φ)*;) where 0 is the zero vector of dimension (N×1) anda(θ, φ) the direction vector associated with the direction (θ, φ), for agiven meshing of the direction of arrival space, estimating the vectorβ_(q,l)(φ_(i)) associated with source i is given by the eigenvector ofthe matrix (A;{tilde over ( )}_(e,q,l)(θ;̂_(i), φ;̂_(i))^(†)A;{tilde over( )}_(e,q,l)(θ;̂_(i), φ;̂_(i)))⁻¹ A;{tilde over ( )}_(e,q,l)(θ;̂_(i),φ;̂_(i))^(†) U;̂_(2q,en)(l) U;̂_(2q,en)(l)^(†) A;{tilde over ()}_(e,q,l)(θ;̂_(i), φ;̂_(i)) that is associated with the minimumeigenvalue.
 8. The method as claimed in claim 7, comprising thefollowing steps: estimating β;̂_(q,l)(φ_(i)), of the vectorβ_(q,l)(φ_(i)) associated with source i by the eigenvector of the matrix(A;{tilde over ( )}_(e,q,l)(θ;̂_(i), φ;̂_(i))^(†)A;{tilde over ()}_(e,q,l)(θ;̂_(i), φ;̂_(i))⁻¹ A;{tilde over ( )}_(e,q,l)(θ;̂_(i),φ;̂_(i))^(†) U;̂_(2q,en)(l) U;̂_(2q,en)(l)^(†) A;{tilde over ()}_(e,q,l)(θ;̂_(i), φ;̂_(i)) that is associated with the minimumeigenvalue, where (θ;̂_(i), φ;̂_(i)) is the estimated direction of sourcei, decomposing the vector β;̂_(q,l)(φ_(i)) into 2^(q−2)(4×1) subvectors,β;̂_(q,l,s)(φ_(i)), 1≦s≦q−2, such thatβ;̂_(q,l)(φ_(i))=[β;̂_(q,l,1)(φ_(i))^(T), . . . ,β;̂_(q,l(q−2))(φ_(i))^(T)]^(T) and ranking the components of each vectorβ;̂_(q,l,s)(φ_(i)) in a (2×2) matrix Γ;̂_(q,l,s)(φ_(i)) such thatΓ;̂_(q,l,s)(φ_(i))[k, j]=β;̂_(q,l,s)(φ_(i))[2(k−1)+j], whereΓ;̂_(q,l,s)(φ_(i))[k, j] and β;̂_(q,l,s)(φ_(i))[k] are the element [k, j]and the component k respectively of the matrix Γ;̂_(q,l,s)(φ_(i)) and ofthe vector β;̂_(q,l,s)(φ_(i)), jointly diagonalizing the matricesΔ;̂_(q,l,s)(φ_(i)), 1≦s≦q−2, to obtain an estimate β;̂(φ_(i)), of thevector β(φ_(i))=;^(Δ)[e^(jφ) _(i) , e^(−jφ) _(i) ]^(T) of the phaseparameters of source i estimate corresponding to the eigenvectorassociated with the maximum eigenvalue, where Δ;̂_(q,l,s)(φ_(i)) isdefined by:Δ;̂_(q,l,s)(φ_(i))=Γ;̂_(q,l,s)(φ_(i))Γ;̂_(q,l,s)(φ_(i))^(†) if q−l=0Δ;̂_(q,l,s)(φ_(i))=Γ;̂_(q,l,s)(φ_(i))Γ;̂_(q,l,s)(φ_(i))^(†) if q−l=1Δ;̂_(q,l,s)(φ_(i))=Γ;̂_(q,l,s)(φ_(i))*Γ;̂_(q,l,s)(φ_(i))^(T) if q−l>1 9.The method as claimed in claim 1, wherein the decomposition stepcomprises a step of estimating the rank of the matrix U;̂_(2q,en)(l). 10.A device for high-resolution direction finding to an arbitrary evenorder, 2q with q>2, in an array comprising N narrowband antennas eachreceiving the contribution from P sources, wherein it comprises aprocessor suitable for implementing the steps of the method as claimedin claim 1.